Contents

Comparisons with CryoSat-2 and PIOMAS sea ice thickness estimates

In this notebook, we compare our ICESat-2-derived sea ice thickness estimates with various estimates derived from ESA’s CryoSat-2 radar altimeter, a combined ICESat-2/CryoSat-2 product and PIOMAS. See the CryoSat-2 wrangling notebook for more insight into the different CryoSat-2 products.

Notes:

  • Doesn’t include spatial comparison stats, but would be pretty trivial to add that in - working on it!

  • ALso lots more to be done in terms of assessing consistency between CS-2 and IS-2 and where/why that breaks down.

Import notebook dependencies

import xarray as xr 
import pandas as pd
import numpy as np
import itertools
import pyproj 
from netCDF4 import Dataset
import scipy.interpolate 
from utils.read_data_utils import read_book_data,read_cs2_book_data  # Helper function for reading the data from the bucket
from utils.plotting_utils import compute_gridcell_winter_means, interactiveArcticMaps, interactive_winter_mean_maps, interactive_winter_comparison_lineplot # Plotting

# Plotting dependencies
import cartopy.crs as ccrs
from textwrap import wrap
import hvplot.xarray
import holoviews as hv
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh # Helps avoid some weird issues with the polar projection 
%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 150 # Sets figure size in the notebook

# Remove warnings to improve display
import warnings 
warnings.filterwarnings('ignore')

Read in the aggregated monthly gridded winter Arctic sea ice data

Below, we’ll read in the aggreated IS-2/CS-2 data and restrict it to the Inner Arctic region.

# Read/download the IS-2 data that includes the various CryoSat-2 products too

book_ds = read_cs2_book_data()

# Restrict to the Inner Arctic domain 
innerArctic = [1,2,3,4,5,6]
book_ds = book_ds.where(book_ds.region_mask.isin(innerArctic))

# Set analysis time range

start_date = "Nov 2018"
end_date = "Apr 2021"
date_range = pd.date_range(start=start_date, end=end_date, freq='MS') # MS indicates a time frequency of start of the month
#date_range = date_range[((date_range.month <5) | (date_range.month > 8))]

Choose CS-2 product

Choose from the following options: AWISMOS, CPOM, GSFC, UBRIS, KK. See the CryoSat-2 wrangling notebook for more insight into the different CryoSat-2 products.

cs2_var='AWISMOS'

Generate a new dataset using a common data coverage mask

Approach is to mask primarily based on the IS2 data then PIOMAS and CS-2 if there is valid data in those datasets fro the given month (really just to stop cases when CS-2 data is completely missing from that month, PIOMAS has complete seasonal coverage).

Likely a much easier way to do this in xarray but couldn’t figure it out |

book_ds_common_mask_list=[]
for date in date_range:
    masked_month = book_ds.sel(time=date).copy(deep=True)
    # Mask everything where missing ICESat-2 data
    masked_month = masked_month.where(masked_month["ice_thickness_int"]>0.0)
    if np.isfinite(np.nanmean(masked_month["cs2_sea_ice_thickness_"+cs2_var])):
        # If there is valid data in the CS-2 dataset then mask all data using this dataset
        masked_month = masked_month.where(masked_month["cs2_sea_ice_thickness_"+cs2_var]>0.0) 
    if np.isfinite(np.nanmean(masked_month["piomas_ice_thickness"])):
        # If there is valid data in the PIOMAS dataset then mask all data using this dataset (I think this never really happens)
        masked_month = masked_month.where(masked_month["piomas_ice_thickness"]>0.0) 
    
    book_ds_common_mask_list.append(masked_month)

book_ds_common_mask = xr.concat(book_ds_common_mask_list, "time")

# Generate monthly means
mean_monthly_data = book_ds.mean(dim=("x","y"), skipna=True, keep_attrs=True)
mean_monthly_data_common_mask = book_ds_common_mask.mean(dim=("x","y"), skipna=True, keep_attrs=True)
#mean_monthly_gsfc = cs2_gsfc.mean(dim=("x","y"), skipna=True, keep_attrs=True)

# Holoviews requires us to generate individual line plots then combine later
is2_lineplot_p = mean_monthly_data.ice_thickness_int.hvplot.line(grid=True, label="IS2SITMOGR4", line_dash='dashed', color='gray', frame_width=700, frame_height=250) * mean_monthly_data.ice_thickness_int.hvplot.scatter(marker='x', color='gray', size=40) # Overlay scatter plot to add markers
is2_lineplot_p_cm = mean_monthly_data_common_mask.ice_thickness_int.hvplot.line(grid=True, label="IS2SITMOGR4 CM",  line_width=1.5,color='k', frame_width=700, frame_height=250) * mean_monthly_data_common_mask.ice_thickness_int.hvplot.scatter(marker='o', color='k', size=40) # Overlay scatter plot to add markers
cs2_lineplot_p = mean_monthly_data['cs2_sea_ice_thickness_'+cs2_var].hvplot.line(grid=True,  line_dash='dashed', label="CryoSat-2", color='g', frame_width=700, frame_height=250) * mean_monthly_data['cs2_sea_ice_thickness_'+cs2_var].hvplot.scatter(marker='x', color='g', size=40) # Overlay scatter plot to add markers
cs2_lineplot_p_cm = mean_monthly_data_common_mask['cs2_sea_ice_thickness_'+cs2_var].hvplot.line(grid=True, label="CryoSat-2 CM",  line_width=1,color='g', frame_width=700, frame_height=250) * mean_monthly_data_common_mask['cs2_sea_ice_thickness_'+cs2_var].hvplot.scatter(marker='o', color='g', size=40) # Overlay scatter plot to add markers
piomas_lineplot_p = mean_monthly_data.piomas_ice_thickness.hvplot.line(grid=True, label="PIOMAS", line_dash='dashed', color='c',  frame_width=700, frame_height=250) * mean_monthly_data.piomas_ice_thickness.hvplot.scatter(marker='x', color='c', size=30) # Overlay scatter plot to add markers
piomas_lineplot_p_cm = mean_monthly_data_common_mask.piomas_ice_thickness.hvplot.line(grid=True, label="PIOMAS CM", color='c', line_width=1,frame_width=700, frame_height=250) * mean_monthly_data_common_mask.piomas_ice_thickness.hvplot.scatter(marker='o', color='c', size=40) # Overlay scatter plot to add markers

r_squared_str = '%.02f' %(xr.corr(mean_monthly_data_common_mask.ice_thickness_int, mean_monthly_data['cs2_sea_ice_thickness_'+cs2_var]).values**2)
mb_str = '%.02f' %(np.mean(mean_monthly_data_common_mask.ice_thickness_int-mean_monthly_data_common_mask['cs2_sea_ice_thickness_'+cs2_var]).values)
sd_str = '%.02f' %(np.std(mean_monthly_data_common_mask.ice_thickness_int-mean_monthly_data_common_mask['cs2_sea_ice_thickness_'+cs2_var]).values)

r_squared_str_pmas = '%.02f' %(xr.corr(mean_monthly_data_common_mask.ice_thickness_int, mean_monthly_data['piomas_ice_thickness']).values**2)
mb_str_pmas = '%.02f' %(np.mean(mean_monthly_data_common_mask.ice_thickness_int-mean_monthly_data_common_mask['piomas_ice_thickness']).values)
sd_str_pmas = '%.02f' %(np.std(mean_monthly_data_common_mask.ice_thickness_int-mean_monthly_data_common_mask['piomas_ice_thickness']).values)


# Combine plots and display
comb_plot = (is2_lineplot_p*is2_lineplot_p_cm*cs2_lineplot_p*cs2_lineplot_p_cm*piomas_lineplot_p*piomas_lineplot_p_cm).opts(hv.opts.Layout(shared_axes=True, merge_tools=True))
comb_plot.opts(title=" IS2SITMOGR4 vs CS-2 ("+cs2_var+r") r²: "+r_squared_str+" "+mb_str+" +/- "+sd_str+r" (m)    IS2SITMOGR4 vs PMAS r²: "+r_squared_str_pmas+"  "+mb_str_pmas+" +/- "+sd_str_pmas+" (m)", fontsize=10)
#comb_plot.opts(annotate="Mean sea ice thickness (meters)")
comb_plot.opts(ylabel="Mean sea ice thickness (meters)")
#comb_plot.opts(annotate=r_str)
comb_plot.opts(xticks=book_ds.time.values[0::6])
display(comb_plot)
hv.save(comb_plot, './figs/icesat-2_cryosat-2'+cs2_var+'.png', fmt='png', dpi=300)

# Interactive map plot for a given month

date='2019-03'

# Use either all data or data on a common mask
data1 = book_ds['cs2_sea_ice_thickness_'+cs2_var].sel(time=date[0:7])
data2 = book_ds.ice_thickness_int.sel(time=date[0:7])
#data1 = book_ds_common_mask['cs2_sea_ice_thickness_'+cs2_var].sel(time=date[0:7])
#data2 = book_ds_common_mask.ice_thickness_int.sel(time=date[0:7])
pl1 = interactiveArcticMaps(data1, vmin=0, vmax=4, colorbar=False, title=date[0:7]+'  CS-2 ('+cs2_var+')')
pl2 = interactiveArcticMaps(data2, vmin=0, vmax=4, colorbar=True, clabel="sea ice thickness (m)", title='ICESat-2') 
diff_pl = interactiveArcticMaps((data2-data1), vmin=-1.5, vmax=1.5, cmap="coolwarm", clabel="difference", title="Difference (b - a)") 
data_pl = (pl1+pl2+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
#data_pl.opts(title="Winter "+str(years[i])+"-"+str(years[i]+1))
display(data_pl)

    

# Interactive map plot for a given month

date='2019-03'

# Use either all data or data on a common mask
data1 = book_ds.piomas_ice_thickness.sel(time=date[0:7])
data2 = book_ds.ice_thickness_int.sel(time=date[0:7])
#data1 = book_ds_common_mask['cs2_sea_ice_thickness_'+cs2_var].sel(time=date[0:7])
#data2 = book_ds_common_mask.ice_thickness_int.sel(time=date[0:7])
pl1 = interactiveArcticMaps(data1, vmin=0, vmax=4, colorbar=False, title=date[0:7]+' PIOMAS')
pl2 = interactiveArcticMaps(data2, vmin=0, vmax=4, colorbar=True, clabel="sea ice thickness (m)", title='ICESat-2') 
diff_pl = interactiveArcticMaps((data2-data1), vmin=-1.5, vmax=1.5, cmap="coolwarm", clabel="difference", title="Difference (b - a)") 
data_pl = (pl1+pl2+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
#data_pl.opts(title="Winter "+str(years[i])+"-"+str(years[i]+1))
display(data_pl)